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1.  Introduction 

Generalized  Reduced  Gradient  (GRG)  Methods  are  algorithms  for  solving 
nonlinear  programs  of  general  structure.  An  earlier  paper  Til  discussed  the 
basic  principles  of  GRG  and  presented  the  preliminary  design  of  a GRG 
computer  code.  This  paper  describes  a modified  version  of  that  initial 
design,  including  the  experiences  that  led  to  the  modifications.  This 
paper  also  is  Intended  to  serve  as  partial  system  documentation.  The 
code  is  compared  computationally  with  an  interior  penalty  function  code, 
and  anticipated  future  work  on  the  algorithm  is  outlined. 
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2.  Brief  Description  of  Generalized  Reduced  Gradient  AlgorJ.thms 

Generalized  Reduced  Gradient  (GRG)  Algorithms  solve  nonlinear  programs 
of  the  form 


minimize 


subject  to  gjCX)  = 0,  1=1,  NEQ 


0 < g^(X)  < UB(N+1),  l=NEQfl,  M 


(1) 


LB(1)  < X^  < UB(1),  1=1,  N 


where  X Is  a vector  of  N variables.  NEQ,  the  number  of  equality  constraints. 


may  be  zero.  The  functions  g^  are  assumed  differentiable. 


There  are  many  possible  GRG  algorithms.  Their  underlying  concepts 
are  described  in  references  . 1_  - 3j.  This  paper  briefly  describes  the 
version  currently  implemented  in  our  code. 

The  user  submits  the  problem  in  the  above  form.  It  is  converted  to 
the  following  equality  form  by  adding  slack  variables  Xjj^j^  ^+M* 


minimize 


8m+1 


subject  to  g^(X)  - Xjj^^  » 0,  1=1,  M 


where 


LB(i)  UB(i),  1-1,  N+M 

LB(i)  = ub(1)  = 0,  i-N+1,  N+NEQ 


(2) 


LB(i)  = 0 i=N+NE(^l,  N+M 


These  last  two  equations  are  the  bounds  for  the  slack  variables.  The  variables 
^1,...,  % will  be  called  "natural"  variables. 


't 


Let  X satisfy  the  constraints  of  (1) , and  assume  that  NB  of  the 
constraints  are  binding  (i.e.  hold  as  equalities)  at  X.  A constraint  g^  is 
taken  as  binding 


if  |g^  - UB(N+i)|  < EPNEWT 

or  |gi  - LB(N+1)|  < EPNEWT 


i.e.  if  it  is  within  EPNEWT  of  one  of  its  bounds.  The  tolerance  EPNEWT  is 

one  of  the  most  critical  parameters  in  the  code.  It  can  be  set  by  the  user, 

-4 

and  has  a default  value  of  10  . 

GRG  uses  the  NB  binding  constraint  equations  to  solve  for  NB  of  the 
natural  variables,  called  the  basic  variables,  in  terms  of  the  remaining 
N-N3  natural  variables  atJ  the  NB  slacks  associated  with  the  binding  constraints. 
These  N variables  are  called  nonbaslc.  Let  y be  the  vector  of  NB  basic 
variables  and  x the  vector  of  N nonbasic  variables,  with  their  values 
corresponding  to  X denoted  by  (y,  x).  Then  the  binding  constraints  can  be 
written. 


g(ypc)  » 0 (3) 

where  g is  the  vector  of  NB  binding  constraint  functions.*  The  basic 
variables  must  be  selected  so  that  the  NB  by  NB  bavsls  matrix 

B = (3gi/9yj) 

is  nonsingular  at  X.  Then  the  binding  constraints  (3)  may  be  solved 
(conceptually  at  least)  for  y in  terms  of  x yielding  a function  y(x),  valid 
for  all  (y,x)  sufficiently  near  (y,x).  This  reduces  the  objective  to  a 


*The  definitions  of  g are  extended  here  to  Include  the  slacks. 
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function  of  x only 

and  reduces  the  original  problem  (at  least  in  the  neighborhood  of  (y,x)), 
to  a simpler  reduced  problem 

minimize  F(x) 

subject  to  Jl  £ X < u 

where  H and  u are  the  bound  vectors  for  x.  The  function  F(x)  is  called  the 
reduced  objective  and  its  gradient,  VF(x),  the  reduced  gradient . 

This  GRG  code  solves  the  original  problem  (1)  by  solving  (perhaps  only 
partially)  a sequence  of  reduced  problems.  The  reduced  problems  are  solved 
by  a gradient  method.  At  a given  iteration  with  nonbasic  variables  x 
and  basic  variables  y,  is  computed,  and  VF(x)  is  evaluated  as  follows: 

u = (9Rj^3^/3y)^B~^ 

9F73xj^  = ~ “ 9g/3xj^ 

A search  direction  d is  formed  from  VF(x)  and  a one  dimensional  search  is 
initiated,  whose  goal  is  to  solve  the  problem 

minimize  F(x  + ad), 
a > 0 

This  minimization  is  done  only  approximately  and  may  be  terminated  for  a 
variety  of  reasons  (see  section  5).  It  is  accomplished  by  choosing  a 
sequence  of  positive  values  {aj^,a2,...}  for  a.  These  are  generated  by 
subroutine  DMINRG,  described  in  section  5.  For  each  value  a^,  F(x+a^d) 
must  be  evaluated. By  (4),  this  is  equal  to  gj^j^(y(x+a^d) , x+a^d) , so  the  basic 
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variables  y(x+a£d)  must  be  determined.  These  satisfy  the  system  of  equations 

g(y,  x4<Xj^d)  » 0 

This  system  is  solved  by  a variant  of  Newtons  method.  If  Newtons  method 
converges,  and  no  constraints  are  violated  at  the  solution,  a new  a value 
is  selected  and  the  one  dimensional  search  process  continues.  If  any 
constraints  or  any  bounds  on  basic  variables  y are  violated, the  code 
determines  a new  a value  such  that  at  least  one  such  new  constraint  or  variable 
is  at  its  bounds  and  all  others  are  satisfied.  If  certain  conditions  are 
met  (see  description  of  subroutine  DMINRG,  section  5) , the  new  constraint  is 
added  to  the  set  of  binding  constraints,  the  one  dimensional  search  is 
terminated,  and  solution  of  a new  reduced  problem  begins. 

3.  System  Overview 

The  GRG  code  described  here  is  composed  of  a main  program  and  a number 
of  subroutines.  It  is  written  in  FORTRAN  IV  and  is  currently  operative  on 
a UNIVAC  1108  at  Case  Western  Reserve  University  and  an  IBM  370-145  at 

Cleveland  State  University,  uses  double  precision  arithmetic. 

System  input  is  described  in  the  user  documentation.  The  only  user 
supplied  subroutine  required  is  GCOMP,  which  ccmpuces  the  functions  g^  for 
given  X.  The  code  requires  first  derivatives  of  the  functions  g^,  but  these 
may  be  computed  by  a system  subroutine  PARSH  using  finite  difference 
approximations.  Alternatively,  the  user  may  supply  a subroutine  PARSH  which 
computes  first  derivatives  by  analytic  formulas  or  other  means. 

The  code  operates  in  two  phases.  If  the  initial  vector  jc  does  not 
satisfy  one  or  more  of  the  g^  constraints,  phase  I finds  an  initial  feasible 
point  or  determines  that  there  is  none.  This  is  done  by  minimizing  a phase  I 


objective  function,  which  is  the  sum  of  the  constraint  violations. 

Phase  II  starts  with  an  initial  feasible  point  and  attempts  to  minimize  the 
aser-supplied  objective  As  with  any  other  NLP  algorithm  any  solution 

found  may  oe  only  a local  rather  than  a global  minimum.  In  phase  I,  this 
means  that  a feasible  point  may  not  be  located  even  if  one  exists.  A popular 
procedure,  if  local  optima  appear  to  be  a problem,  is  to  try  a variety  of 
starting  points.  If  the  same  final  point  is  obtained,  it  is  likely 
that  this  is  a global  solution.  A suggested  algorithm  for  generating 
alternative  starting  points  is  described  in  reference  . 

y 

4 . Subroutines  Comprising  the  Code 

The  current  GRG  code  is  composed  of  a main  program,  MAINRG,  and  11 
subroutines.  These  are  described  briefly  in  the  following  subroutine 


dictionary. 


Subroutine  Dictionary 


Subroutine 

Name 

1.  MAINRG 

2 . SUMRY 

3.  GCOMP 

4.  PARSH 

5 . GRG 

6.  DMINRG 

7.  REDOB J 


Function 

Calls 

Called 

by 

Not  a subroutine.  Reads,  edits. 

GRG 

and  prints  input. 

SUMRY 

— — 

Optional  user  supplied  subroutine 
which  prints  additional  solution 
output,  beyond  that  provided  by 
GRG 

MAINRG 

User  supplied  subroutine.  Given 

— 

GRG 

current  X vector,  computes  vector 
of  ^H-l  function  values  G,  where 
G(l),...,  G(M)  are  constraint 
function  values  and  G(M+1)  is  the 
objective 

NEWTON 

Given  current  X and  G vectors. 

If  not 

computes  array  GRAD(Iifl,  N), 

user 

CONSEf 

whose  (I,J)  element  is  the 

supplied. 

partial  derivative  of  G(I)  with 

calls 

respect  to  X(J).  May  be  user 
supplied.  If  not,  there  is  a 
system  subroutine  PARSH  which 
computes  GRAD  by  foi'ward 
difference  approximation 

GCOMP 

Controls  main  iterative  loop. 

GCOMP 

MAINRG 

Computes  initial  BINV  and 

CONSBS 

search  direction.  Calls  one 

REDGRA 

dimensional  search  subroutine 
DMINRG.  Tests  tor  optimality, 
updates  H matrix 

DMINRG 

Performs  one  dimensional 

. REDOBJ 

GRG 

search 

CONSBS 

Computes  values  of  basic 

NEWTON 

DMINRG 

variables  for  given  values  of 

GCOMP 

nonbasics  by  calling  NEWTON. 

CONSBS 

Takes  action  if  NEWTON  doesn't 

PHIOBJ 

converge.  Checks  for  constraint 
violations.  If  any  are  violated, 
finds  feasible  point  where  some 
Initially  violated  constraint  is 
binding  and  others  satisfied. 


P!PP9PP!9 
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Subroutine  Dictionary  - p.2 


Subroutine 

Name 


Function 


Calls  Called 

by 


8.  NEWTON  Uses  Newtons  Method  to  compute 

values  of  basic  variables  for 
given  values  of  nonbasics.  If 
convergence  not  achieved,  sets 
flag  and  returns 

9.  CONSBS  CotujHites  Basis  Inverse,  BINV 


UCOMP  REDOBJ 


PARSH 


GRG 

REDOBJ 

DMINRG 


10.  REDGRA 


11.  PUIOBJ 


12.  TAI^r. 


Given  BINV  and  GRAD,  computes 
Lagrange  multiplier  vector  U, 
and  reduced  gradient  of  either 
phase  I or  phase  II  objectives, 
GRADF 

If  any  constraints  are  violated, 
computes  phase  I objective, 
equal  to  the  sum  of  constraint 
violations.  Stores  this  as 
G(M+1),  stores  original  G(M+1) 
as  TRU03J. 

Computes  tangent  vector  v as  v = 


GRG 


REDOBJ 


GRG 


(8g/9x)d 


5.  Subroutine  Flow  Charts  and  Descriptions 

The  flow  charts  In  this  section  are  In  aggregated  rather  than  detailed 
form.  Their  purpose  Is  to  describe  overall  program  logic.  However,  they 
correspond  fairly  closely  to  the  acrual  FORTRAN  code.  In  particular,  all 
array,  variable,  and  subroutine  names  used  aiethe  same  as  In  the  code. 
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Subroutine  GRG 

The  subroutine  begins  by  calling  CONSBS  to  invert  the  initial  basis.  If 
NCAND  is  not  zero  in  block  1,  the  user  has  specified  an  initial  candidate 
list  for  CONSBS.  Otherwise,  all  variables  are  candidates.  The  arrays 
lABOVE  and  IBULOW  in  block  2 are  the  sets  of  Indices  of  constraints  which 
vie i ate  their  upper  and  lower  bounds  respectively. 

The  Broyden  - Fletcher  - Shanno  (BFS)  variable  metric  algorithm  181  is 
used  to  generate  the  search  direction  d.  This  method  uses  an  N x N matrix,  H. 

In  block  3,  H is  Ir  tlalized  to  a 4iagonal  matrix  with  diagonal  element  zero 
if  the  1 nonbaslc  variable  is  at  a bound  and  unity  otherwise.  The  test  in 
block  4 is  true  if  either  of  two  optimality  tests  are  passed.  The  first  test 
checks  if  the  following  conditions  are  met; 

for  i = 1,N  but  x^  not  a slack  variable 

for  an  equality  constraint 

x(l)  = LB  (1)  ^ GRADF(i)  > -EPSTOP 

x(i)  = UB  (i)  ^ GRADF(i)  < EPSTOP 
LB(i)  < x(i)  <UB(i)r^|GRADF(i)|<  EPSTOP 

The  quantities  x^  are  the  current  nonbasic  variables  and  GRADF(i)  is  the 
til 

i component  of  the  reduced  gradient,  (see  section  2).  This  tests  whether  the 
Kuhn-Tucker  optimality  conditions  [7]  are  satisfied  to  within  EPSTOP,  a small  positive 
number  which  can  be  controlled  by  the  user.  The  slack  variables  for  equality 
constraints  (l.e.  the  variables  X(N+1)  to  X(NvNEQ))  are  excluded  from  the  test 
because  they  must  be  zero  in  any  feasible  solution. 


I 

I 
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The  second  optimality  test  checks  if  the  condition 
ABS(FM  - OBJTST)  < EPSTOP  * ABS(OBJTST) 


I 


1 

-itn 


i . 


is  satisfied  for  NSTOP  consecutive  iterations.  In  the  above,  FM  is  the 

current  objective  value  and  OBJTST  is  the  Objective  value  at  the  start  of  the  previous 

one  dimensional  search. 

There  are  two  tests  for  resetting  H in  block  5.  The  first  tests  whether 
the  scalar  product 


Z d(i)  * GRADF(i) 
i 

is  negative.  If  not,  the  search  direction  d(i)  will  not  yield  an  immediate 
decrease  in  the  objectiv*,  and  H must  be  reset.  This  condition  can  occur 
due  to  numerical  error  in  computing  H or  to  Inaccuracies  in  the  one  dimensional 
search.  The  second  test  checks  if 

max  |d(i)  | < 10  ^ 

i.e  if  d is  too  small.  Neither  of  the  latter  two  tests  can  be  true  immediately 
after  H is  reset. 

In  block  6,  AI.FMAX  is  the  largest  value  of  a for  which  each  component  of 
X + ad  satisfies  its  bounds.  The  tangent  vector,  V,  in  block  7 is  used  to 
compute  initial  values  for  the  basic  variables  in  subroutine  REDOBJ. 

The  variable  IFLAG  in  block  8 is  set  to  3 either  in  REDOBJ,  if  phase  I 
ends,  or  in  DMINRG,  if  a new  binding  constraint  is  to  be  added  to  the  basis. 
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In  either  case  a new  reduced  problem  Is  to  be  solved,  so  H is  reinitialized 
and  the  entire  procedure  begins  again.  IFLA6  is  set  to  6 either  in 
DMINR6,  if  too  many  NEVTTON  iterations  have  been  taken,  or  in  REDOBJ, 
if  NEWTON  falls  to  converge.  A new  one  dimensional  search  is  initiated 
but  the  reduced  problem  remains  the  same,  so  H is  not  reset. 

In  block  10,  page.  2,  the  binding  constraints  are  checked  to  see  if  any  have 
become  strictly  satisfied  during  the  one  dimensional  search.  If  so,  a new, 
smaller  basis  Inverse  is  constructed  in  CONSBS,  and  a new  reduced  problem 
solution  begins. 

If  the  one  dimensional  search  ends  with  IFL'AG  = 0,  corresponding  to 
an  unconstrained  minimum  being  located  along  the  search  direction  d,  H is 
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I 

i 

I 
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updated  and  a new  iteration  begins.  The  update  used  depends  on  whether  or 
not  any  nonbasic  variable  has  reached  a bound,  i.e.  if  the  step  size,  ALPH, 
is  equal  to  ALFMAX.  The  updating  formulas  in  blocks  11  and  12  are  given 
in  reference  1^5]  . Since  H is  a symmetric  matrix,  only  the  diagonal  and 
super-diagonal  elements  are  needed,  and  these  are  stored  in  a linear  array 
H in  row  order. 


& 


15 


Subroutine  I»1INRG 

Subroutine  GRG  provides  search  directions  for  the  one  dimensional  search 
subroutine,  DMINRG, in  which  the  variables  of  the  problem  are  assigned  new  values. 

This  subroutine  finds  a first  local  minimum  for  the  problem 
minimize  F(x  + ad) 

The  direction  d is  always  a direction  of  descent,  i.e. 
d^7F  (5)  < 0 

This  subroutine  searches  for  three  a values.  A,  B,  and  C,  which 
^atlsfy 

0 _<  A < B < C 
F (x  + Ad)  > F (x  + Bd) 

and 

F (x  + Cd)  > F (x  + Bd) 

Then  the  interval  [A,C]  contains  a local  minimum  of  F (x  + ad).  In  block  11 
of  the  flow  chart,  a quadratic  in  a is  passed  through  A,  B,  and  C,  with  itc 
minimum  at  D.  The  point  D is  taken  as  an  estimate  of  the  optimal  a and  z return 
is  made. 

In  finding  (A,B,C)  the  choice  of  initial  step  size,  a^,  (block  1,  page  1), 

is  important.  With  Goldfarbs  algorithm  or  other  variable  metric  methods,  a^ 

is  set  equal  to  the  optimal  a value  from  the  previous  search  except  when  this 

causes  too  large  a change  in  the  variables . The  theoretical  basis  for  this  is  that, 

as  a variable  metric  converges,  the  optimal  a values  should  converge  to  1,  the 

optimal  step  size  for  Newton's  Method.  Hence  the  previous  optimal  step  is  a good 

approximation  to  the  current  one.  This  must  be  modified  when  the  method  is  restarted, 

for  example  when  a new  constraint  is  encountered  or  the  basis  is  changed,  since  then 

an  optimal  step  much  less  than  unity  is  generally  taken.  Hence,  we  require  that  the 

~3 

change  in  any  nonbasic  variable  larger  than  10  in  absolute  value  not  exceed  .05 
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times  its  value,  while  the  change  in  any  variable  smaller  than  10  in  absolute 
value  cannot  exceed  0.1.  If  the  largest  o value  meeting  these  conditions  is  a^, 
and  is  the  step  size  found  by  DMINRG  at  interation  i,  then  is  equal  to 
a = min  (a . , . a^) 

if  the  previous  search  terminated  with  an  interpolation,  and  = a'*'  otherwise. 

The  loop  2 - 3 - 4 - 5 halves  the  step  size  until  a value  FB  < FA  is  achieved, 

or  until  LOOPCT  ■ 10.  The  variable  IFLAG  in  block  3 is  set  in  REDOBJ  - to  3 if  a 

new  constraint  or  bound  on  basic  variable  was  encountered  and  a new  basis  was 

constructed,  and  to  6 if  either  NEWTON  call  in  REDOBJ  did  not  converge. 

The  test  in  block  6 of  the  one  dimensional  search  flow  chart  is  false  only  if 

the  step  size  has  been  halved  at  least  once  in  2-  3-4-  5,  in  which  case  Ki  is  the 

function  -;alue  corresponding  to  C.  It  also  insures  that  the  subroutine  will  cut 

back  the  subroutine  will  cut  back  the  step  size  if  a large  function  value  is  returned 

by  REDOBJ.  This  is  used  to  force  a cutback  when  the  Newton  algorithm  in  REDOBJ  does 

30 

not  converge  and  an  improved  point  has  not  been  found,  by  setting  FB  to  10 

The  loop  7-8-9-10  doubles  the  step  size  each  time  until  the  points  A,  B, 

C bracket  the  minimum.  The  test  in  block  7 is  true  if  the  NEWTON  algorithm  in  REDOBJ 
took  more  than  5 iterations  to  converge.  Experience  has  shown  that,  in  this  case, 
the  next  step  with  C 2B  is  almost  certain  not  to  converge  in  the  limit  of  10 
iterations.  This  test,  and  related  logic,  has  reduced  overall  computational  effort 
significantly  - see  sections  6 and  7. 

Subroutine  DMINRG  also  includes  logic  to  insure  that  the  step  taken  is  no  larger 
than  ALFMAX.  To  simplify  the  exposition,  this  logic  has  not  been  included  in  the 
flow  chart.  Before  returning,  DMINRG  picks  up  the  best  objective  value  encountered 
during  the  search.  This  is  done  in  block  12,  page  2.  The  quantities  XBEST,  GBEST, 
ALFBST  are  computed  in  subroutine  REDOBJ,  which  will  now  be  described. 
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Subroutine  REDOBJ  - Introduction 

This  subroutine  (REDuced  OBJectlve)  is  equivalent  to  the  subroutine 
which  evaluates  the  objective  function  in  a procedure  for  unconstrained 
minimization.  However,  it  is  much  more  complicated.  REDOBJ  is  called  from 
the  one  dimensional  search  subroutine  DMINRG.  Prior  to  the  call  to  REDOBJ, 
DMINRG  chooses  a value  for  the  step  size  a,  and  computes  nev;  values  of  the  non- 
basic  variables,  x,  equal  to  x + ad  (d  is  the  search  direction).  Then 
REDOBJ  is  called.  It  attempts  to  compute  the  corresponding  values  of  the 
basic  variables,  by  solving  the  system  of  NB  nonlinear  equations 
8i  <.y»  + ad)  = 0 lelBC 

for  the  NB  basic  variables  y,  vdiere  IBC  is  the  index  set  of  binding 
constraints.  This  system  is  solved  by  subroutine  NEWTON. 

If  a solution  is  obtained  then  all  of  the  constraints  are  checked  to 
see  that  none  are  violated.  If  not,  and  if  no  new  constraints  are  binding 

then  the  current  objective  value  is  compared  to  the  previous  best  value  • 

If  it  is  lower,  the  current  values  of  the  variables,  constraints, 

objective  and  step  size  are  stored  as  XBEST,  GBEST  and  ALFBST 
before  returning  to  the  calling  subroutine  (DMINRG). 

If  the  pseudo-Newton  algorithm  (NEWTON),  used  to  solve  for  the  basic 
variables  in  terms  of  the  nonbasics,  does  not  converge  then  one  of  two 
alternatives  is  taken.  If  at  least  one  improved  point  had  been  fcund  during 
the  linear  search  process  in  DMINRG,  then  the  best  such  value  found  is 

accepted  as  the  optimum  in  this  search  direction  and  the  search  terminated. 

If  no  improved  point  had  been  found  (i.e.  ALFBST  = 0)  then  the  objective 
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function  is  assigned  a large  value  (G(M  + 1)  = 10  ) to  force  the  linear  search 

process  to  cut  back  the  step  size  ALPH  and  to  try  again. 

If,  after  the  NEWTON  process  has  converged,  one  of  the  constraints  is 
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violated  then  REDOBJ  attempts  to  find  the  largest  value  of  ALPH  such  that 
no  coi.straints  are  violated  and  at  least  one  new  constraint  is  binding. 

If  successful,  control  is  returned  to  DMINRG  where  a new  search  is  initiated. 
If  not  successful  then  the  same  action  is  taken  as  if  NEWTON  did  not  converge. 


f 
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REDOBJ  - Detailed  Description 

Given  the  new  values  for  the  nonbasic  variables  x in  terras  of  their  previous 
values  X,  the  current  search  direction  d and  a step  size  a then  the  new  values  for  the 
basic  variables  y are  determined  by  solving  the  system  of  non-linear  equations 
gj.  (y»  X + ad)  = 0 ielBC 

where  IBC  is  the  index  set  of  binding  constraints.  As  in  [2]  and  [3]  this  is 
accomplished,  in  subroutine  NEWTON,  using  the  pseudo-Newton  algorithm 

y^+i  = y^.  " ^ 

where  g^  is  the  vector  of  binding  constraints.  The  algorithm  is  called  pseudo-Newton 
because  B is  evaluated  once  at  the  initial  point  of  the  search,  X,  instead  of  being  re- 
evaluated at  each  step  of  the  algorithm,  as  in  the  standard  Newton  method. 

An  initial  estimate  of  the  solution  is  computed  by  linear  extrapolation  in 
block  1 on  page  1 of  the  REDOBJ  flow  chart.  Consider  the  tangent  plane  to  the 
constraint  surface  at  X.  This  is  the  set  of  all  vectors  (a,b)  satisfying 
(3g/3y)a  + (3g/3x)  b = 0 

where  all  partial  derivative  matrices  are  evaluated  at  X.  In  GRG,  the  change  in 
X,  b,  is  given  by 
b = ad 

The  corresponding  vector  a is  called  the  tangent  vector,  V.  Since  any  scale 
factor  multiplying  V is  unimportant,  we  may  as  well  take  a =*  1,  yielding 
V = - (3g/3y)"^  (3g/3x)d 

In  our  program,  V i.,  computed  at  X,  the  initial  point  of  the  one  dimensional  search 
in  block  7,  pg.  2,  of  subroutine  GRG.  In  block  1 of  REDOBJ,  this  vector  is  used 
to  find  initial  values,  y^,  by  the  formula 
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Yo  = y + V 

Using  these  initial  values,  Newton  finds  the  feasible  point  Then,  at 

V is  not  recomputed.  The  old  V is  used,  but  emanating  now  from  Xj^,  to 
yield  the  next  set  of  initial  values  as 
Yo  = Yl  + ■ °‘l^^ 

Using  these,  Newton  finds  a new  point  X2*  This  procedure  is  repeated  until 

Newton  falls  to  converge  or  until  the  one  dimensional  search  is  over. 

Newton  is  considered  to  have  converged  i 2 the  condition 

NORMG  = max  |g.(X^.)|  < EPNEWT 
lelBC 

is  met  within  ITLIM  interatlons.  Currently  EPNEWT  = 10“^  and  ITLIM  = 10. 

If  NORMG  has  not  decreased  from  its  previous  value  (or  the  above  condition 
is  not  met  in  10  Iterations)  Newton  has  not  converged.  In  block  2,  ALFBST 
is  the  value  of  a at  the  best  point  found  thus  far  in  the  one-dimensional 
search.  If  ALFBST  = 0,  i.e.  no  better  point  has  been  found,  the  objective 
is  set  to  10^®  in  block  3.  This  forces  DMINRG  to  cut  back  the  step  size. 

The  Newton  algorithm  will  converge  if  the  step  size  is  "small  enough".  If 
ALFBST  > 0,  at  least  one  improved  objective  value  has  been  found  by  DMINRG, 
so  CONSBS  is  called  and  the  search  is  terminated  by  setting  IFLAG  to  6 in 
block  4. 

Once  Newton  has  converged,  we  check  for  constraint  violations  on  page  2. 
There  are  a number  of  reasons  why  the  current  step  a may  be  too  large; 

(1)  A strictly  satisfied  constraint  (one  with  index  in  the  set  IRC) 
mav  have  violated  an  upper  or  lower  bound. 

(2)  A constraint  in  lABOVE,  the  set  of  constraints  initially  violating 
their  upper  bounds,  may  violate  a lower  bound. 

(3)  A constraint  in  IBELOW  may  violate  an  upper  bound. 
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(4)  A basic  variable  may  violate  a lower  or  upper  bound.  j 

J 

s 

If  one  or  more  of  these  cases  hold,  a is  reduced  lo  a value,  ALFSTR,  where  | 

t 

no  constraints  are  violated  and  at  least  ona  new  constraint  is  equal  to  a ■ 

bound.  To  determine  this  constraint,  an  estimate  is  made  of  ALFSTR  in  block 

6 for  cases  1,  2,  and  3.  Linear  interpolation  between  the  current  and 

previous  values  of  the  violated  constraint  is  used.  Case  4 is  dealt  with 

in  the  same  way  in  block  8.  If  none  of  the  4 cases  holds  (Y  branch,  block  9), 

the  test  in  block  10  checks  if  we  are  still  in  phase  I (NABOVE  is  the  number 

of  indices  in  lABOVE,  the  set  of  constraints  violating  upper  bounds,  and 

similarly  for  NBELOW) . If  still  in  phase  I (N  branch,  block  10)  block  11 

checks  to  see  if  all  violated  constraints  are  satisfied.  If  so,  phase  I ends 

and  a return  is  made. 

I any  of  the  4 cases  holds,  the  logxc  at  the  top  of  page  3 decides 
whether  case  4 or  one  of  cases  1 - 3 is  to  be  dealt  with.  Assuming  cases 
1 - 3 as  an  example,  we  then  wish  to  solve  the  system 
g;{^  (y(ct),  X + ad)  = 0,  ielBC 
El  (yCa) , X + ad)  = 0 
The  Jacobian  for  this  system  is 

I ' * 

I B 'cl 

- -J 

D4  w 

I 


where 

D4  = Sg^/Sy 
w = (9gj^/3x)'^d 

and  c is  an  NB  component  column  vector  whose  elements  are  (9g^/9x)'^d  for 


ieIBC.  Since  J involves  only  adding  a border  to  the  current  basis  macrix, 


B,  its  inverse  is  easily  computed  if  is  knovm.  This  is  done  in 
block  13.  The  call  to  NEWTON  in  b3ock  14  then  attempts  to  solve  the  system. 

If  NEWTON  falls  to  converge,  the  same  logic  as  described  previously  if; 
applied.  If  it  does  converge,  we  return  to  page  2 to  see  if  any  of  cases 
1-4  still  holds. 

This  procedure  for  satisfying  violated  constraints  terminates  with 
the  Y branch  in  block  9.  If  a basic  variable  has  been  set  to  one  of  its 
bounds,  its  index  is  stored  as  LVl  in  block  12,  p.  3.  In  block  15,  page  4 LV  Is 
replaced  by  LVl.  If  both  LV  and  NNEW  are  zero,  the  set  of  strictly 
satisfied  constraints  is  checked  to  see  if  any  constraints  in  it  are  binding. 
Then  the  best  point  is  updated  if  necessary  and  a return  is  made. 
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Subroutine  CONSBS 

This  subroutine  selects  a set  ot  basic  variables  and  computes  the  basis 
Inverse,  BINV.  Its  Input  Is  a list  of  Indices  of  variables,  the  candidate 
list  ICAND.  The  outputs  of  CONSBS  are  (a)  a new  list  of  binding  and 
strictly  satisfied  constraint  Indices  (b)  a new  list  of  basic  variable 
Indices  and  (c)  the  new  basis  Inverse,  BINV.  In  block  1 of  the  flow  chart 
the  array  IR£M  contains  the  list  of  rows  of  TAB  which  remain  to  be 
pivoted  In.  The  subroutine  operates  In  2 modes.  Indicated  by  the  variable 
MODE.  When  MODE  « 1,  CONSBS  will  choose  pivot  columns  from  whatever 
candidate  list  was  Input  to  It.  If  a basis  Inverse  could  not  be  constructed 
from  columns  In  this  candidate  list,  or  If  the  original  candidate  list 
Included  all  variables  (NCAND  « N,  block  2),  MODE  Is  set  to  2,  and  CONSBS 
will  choose  pivot  columns  from  the  list  of  all  admissible  columns.  A 
column  is  admissible  if  it  is  not  scheduled  to  leave  the  basis  and  if  it  has 
not  yet  been  pivoted  in. 

The  main  loop  of  CONSBS  begins  at  block  3.  A pivot  row  is  chosen  as 
in  block  4.  The  choice  of  ISV  in  block  5 is  motivated  by  the  desire 
co  have  basic  variables  as  far  from  their  bounds  as  possible,  so  that  fewer 
basis  changes  will  be  required.  The  other  criterion  Influencing  the  choice 
of  basic  variables  Is  that  the  basis  matrix  should  be  well-conditioned. 

We  try  to  Insure  this  by  choosing  as  a prospective  pivot  column  that  Index, 
I,  In  ISV  yielding 

max  I TAB  (IROW,  I)| 

lelSV 

This  is  done  in  block  6.  If  the  element  chosen  passes  2 tests  we  nivot 
on  it  (block  8) , transforming  the  Jacobian  and  entering  that  column  into 
BINV.  The  index  of  the  column  pivoted  in  is  stored  in  the  basic  variable 
list,  IBV  (block  7),  and  the 


procedure  is  repeated  for  each  binding  constraint  until  either  BINV 
has  been  constructed  (N  branch,  block  9)  or  the  candidate  list  has  been 
eidiausted  (Y  branch,  block  10^  page  2) . 

The  two  tests  that  a pivot  element  must  pass  are  (a)  Its  absolute 
value  must  be  larger  than  EPSPIV  (currently  lO"^)  and  (b)  the  absolute 
value  of  the  ratio  of  all  other  elemen:s  In  the  pivot  column  to  the 
pivot  element  must  be  less  than  RTOL  (currently  100).  The  first  test 
Insures  that  we  do  not  pivot  on  an  element  that  is  essentially  zero,  while 
the  second  protects  against  the  generation  of  elements  of  large  absolute 
value  in  BINV.  Such  values  are  symptomtic  of  ill-conditioned  basis  matrices. 

It  either  test  is  failed  while  MODE  ■ 1,  we  simply  do  not  pivot  in  the 
current  row,  and  move  on  to  the  next  one. 

If,  when  mode  1 terminates,  BINV  has  not  yet  been  constructed,  we 
attempt  to  complete  its  construction  by  considering  columns  not  in  the 
originalcandidate  list.  In  block  11,  page  2,  the  candidate  list,  ICAND,  is 
reset  to  the  set  of  all  remaining  admissible  columns,  MODE  is  set  to  2,  and 
we  return  to  the  start  of  the  main  iterative  loop.  If,  in  this  second  phase, 
a pivot  element  fails  the  absolute  value  test,  we  temporarily  mark  all  columns 
in  ISV  inadmissible  (by  setting  their  indicator  in  the  IGNORE  array  to  1,  in 
block  12,  page  3)  and  choose:  a new  ISV  array.  If,  in  mode  2,  a column  fails  the 
RTOL  test  (block  14,  page  3),  its  IGNORE  indicator  is  set  to  one.  If  all  admisslbl 
matrix  elements  in  a row  fail  one  of  these  two  tests  (N  branch,  block  13,  page  2) 
we  are  forced  to  consider  variables  within  EPBOUN  of  their  bounds.  This 
is  done  in  block  15,  p.4.  If  all  these  fail  the  tests,  the  slack  variable 
corresponding  to  row  IROW  is  entered  into  the  basis. 
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Subroutine  CONSBS  - p«  1 


r 


START 

_“t 


# 


Initialize;  ICT  0,  IGNORE  (I)  - 0,  I - 1,  N 

' 'ft'.,  :: 


Determine  indices  of  binding  and 
strictly  satisfied  constraints 
NB  «=  number  of  binding  constraints 


i. 

[call  par^ 

TT“"- 

-•  NB  » 0 


fts  1 p.3 


Store  gradients  of  binding  constraints  in  array  TAB 

ft  ■ 


Sort  variables  in  order  of  increasing  Z(J)  where 
Z(J)  - min  {(X(J)  - LB  (J)),  (UB(J)  - X(J))} 


.-ftD 


V 


IREM  (1,2,...,NB) 
NREM  NB 


i' 


J 

1 


LV  » 0 

YN 


N 

V 


Label  variable  leaving  basis 
inadmissible 


V 


-<• 


I 


i 1 5 page 
2 
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'1 
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Subroutine  CONSBS  - p.2 


j 


NCAND  « N 



N 


MODE  2 

t 


1 


start  of 
main  loop 


r~V 


I 




-j 


MODE  *»-  1 


(J> 

ICT  ICT  + 1 


MODE  = 2 

YY 


FOR  I = 1,N,  IGNORE  (I)  -*•  0 


(lo} [_ 


(10 ./ 

ICT  > NREM 


<■ 


N 


...  © 

I IROW  - IREM  (ICT) 




/K 


[ ICT  ■ 
I MODE 


0 

2 




NREM  number  of  rows  not  yet 
pivoted  in 

IREM  l?3t  of  rows  not  yet 
pivoted  in 

ICAND-*-  list  of  admissible 
variables 

NCAND size  of  ICAND 


© 


Select  (up  to)  5 admissible  candidate  variables  with  IGNORE  (I)  = 0 
which  have  Z(J) > EPBOUN  and  with  largest  values  of  Z(J).  Store 
indices  in  ISV.  Let  NSV  “ number  cf  indices  in  ISV 


Scnn  row  IROW  of  TAB.  From  columns  with  indices  in  ISV, 
pick  one  with  element  of  largest  absolute  value.  Let 
element  value  ■ PIV,  column  index  * ICOL 
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Siibrotttlne  CONSBS  ■>  p.  3 


|PIV|  > EPSPIV 


MODE  - 1 


&>'■ 


j For  each  I in  ISV  set  IGNORE  (I)  - 1 > 

.-’i'-i 

1^.2 

1 

For  I - 1,  HB  but  I ^ IROW,  R(I)  ■ -TAB(I,ICOI.)/PIV 

;-  ^ 

toy  R(I)  > RTOl  >4~y^DE  - 1 ^-^""[T* 

- <3>  h 

Store  ICOL  in  list  of  basic  variables^  IBV,  in  - IGNORE  (ICOL)  ■ 1 | 

position  IROW.  Replace  column  ICOL  of  TAB  by  “ — •* 

...unit  r*n 

^ J P-2 


Pivot  on  TAB  (IROW,  ICOL)*  Update  all  columns  of  TAB 


1 Any  rows  remaining  to  pivot  in 

-> v" 


Reorder  coluranTof  T^' v/ith  TiidiTer'in  ba8rc"'variabTrii7t7~ 
IBV,  Into  first  NB  columns  of  TAB.  Order  is  as  specified  in 
IBV.  Since  TAB  and  BINV  occupy  same  storage  locations  (by 
EQUIVALENCE  statement)  these  columns  comprise  BINV 

Construct  index  set  of  nonbasic  variables 


1 

1 


RETURN 
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6*  Changes  in  the  Algorithm 

The  subroutines  of  section  5 have  been  changed  significantly  from  their 
orginal  versions,  described  in  reference  [1].  Here  we  discuss  some  of  these  changes 
and  give  reasons  for  them. 

Subroutine  PARSH 

A finite  difference  version  of  PARSH  has  been  added.  This  computes  partial 

derivatives  of  g^,..,g  , by  simple  forward  differencing,  using  a constant  increment 

X jj  + X 

-4 

of  10  in  each  variable.  It  frees  the  user  of  having  to  code  his  own  PARSH,  a task 
which  can  require  man  - months  of  effort  for  some  problems.  In  comparing  solution 
by  using  analytic  and  finite  difference  derivatives,  little  or  no  degradation  in 
accuracy  or  speed  has  been  noted. 

The  form  of  the  user  ~ supplied  PARSH  has  also  been  changed.  Previously,  the 
code  assumed  that  only  the  gradients  of  currently  binding  constraints  would  be  computed 
by  PARSH.  Now,  gradients  of  all  constraints  are  computed.  This  significantly 
simplifies  the  preparation  of  PARSH.  The  older  version  required  that  PARSH  compute 
the  gradients  of  an  r.rbitrary  subset  of  constraints,  which  is  more  complex  to  code 
than  simply  computing  all  of  them.  The  rest  of  the  program  is  also  simplified.  No 
distinction  need  now  be  made  between  finite  difference  and  analytic  derivatives; 
previously,  finite  difference  would  compute  all  gradients,  analytic  only  those  for 
binding  constraints. 

Subroutine  GRG 

The  Broyden  - Fletcher  - Shanno  (BPS)  variable  metric  algorithm  [8]  has  repi«ced 
the  Davidon  - Fletcher  - Powell  [6]  method.  This  was  done  because  recent  computational 
experience  indicates  that  BFS  is  the  best  of  the  variable  metric  algorithms.  Reference 
[9]  shows  that  BFS  is  less  sensitive  to  errors  in  the  1 dimensional  search  than  DFP, 
while  [10]  provides  evidence  that  periodic  restarting  of  BFS  is  undesirable.  Hence 
a simpler  1 dimensional  search  is  now  used,  and  the  test  in  block  5 of  GRG  does  not 
include  resetting  H periodically.  ' 
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Subroutine  DMINRG 

The  cubic  interpolation  section  has  been  deleted.  This  followed  the  quadratic 
interpolation  (block  11),  and  fit  successive  cublcs  through  4 points  until  certain 
stop  criteria  were  met.  It  was  removed  because: 

(1)  In  about  10  test  problems,  it  was  found  that  deleting  the  cubic  Interpolatlo 

increased  solution  effort  little.  If  at  all,  and  often  decreased  it.  There  are  a 

number  of  reasons  for  this.  Foremost  among  them  is  the  fact  that,  to  compute  the 

reduced  objective  F(x)  exactly,  the  basic  variables  y(x)  must  be  determined  exactly, 

i.e.  each  binding  constraint  must  be  exactly  equal  to  zero.  Of  course  this  is 

impossible  in  practice,  and  the  binding  constraints  are  only  within  EPNEWT  of  their 

-4 

bounds,  with  the  default  value  of  EPNEWT  currently  equal  to  10  . Hence  F(x)  is 

computed  significantly  less  accurately  than  if  the  problem  were  unconstrained,  perhaps 
only  to  4 or  5 significant  figures.  With  such  noise  in  the  function  evaluation,  cubic 
interpolation  rarely  achieves  Improved  function  values. 

(2)  There  is  much  evidence  in  unconstrained  ttinimizatlon  that  a "sloppy"  one 
dimensional  search  achieves  better  overall  results  than  a more  exact  one. 

(3)  The  code  was  significantly  simplified  and  shortened. 

The  current  interpolation  strategy  is  simply  to  bracket  the  minimum,  fit  a 
quadratic,  choose  the  best  of  all  points  evaluated  during  the  current  search,  and 
return.  Multiple  interpolations  are  no  longer  performed.  Choosing  the  best  point 
is  necessary,  because  the  D point  selected  by  the  quadratic  interpolation  (block  11, 
DMINRG  flow  chart)  may  not  be  the  best,  and  choosing  it  can  yield  a final  objective 
value  worse  than  the  value  at  the  start  of  the  search. 

Two  new  checks  have  been  added  to  the  quadratic  interpolation.  The  first  test 
requires  that  the  product  (A-B)  * (B-C)  * (C-A)  exceed,  in  absolute  value,  a tolerance 
EPSQl,  while  the  second  requires  that  the  three  points  (A, FA),  (B,FB)  and  (C,FC) 
enclose  an  area  no  less  than  a tolerence  EPSQ2.  These  tests  ensure  that  the  quadratic 
interpolation  is  numerically  stable,  and  no  interpolation  is  performed  if  either  of  thi 


fails. 
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Subroutine  REDOBJ 

The  actions  taken  when  NEWTON  does  not  converge  (see  pages  1 and  3 of  the  flow 

chart)  are  new.  In  the  current  logic,  if  NEWTON  does  not  converge  then  either  the 

best  point  found  so  far  is  accepted  as  the  minimum  and  the  one  dimensional  search 

terminated  or,  if  no  such  point  has  been  found,  the  objective  value  at  the  current 
30 

point  is  set  to  10  in  order  to  force  the  step  size  to  be  cut  back  in  DMINRG. 

Previous  logic  tried  to  force  convergence  by  first  computing  an  approximate  B ^ and, 
if  that  also  failed,  then  an  exact  B both  evaluated  at  the  last  feasible  point. 

The  approximate  B ^ was  a dismal  failure,  and  even  the  exact  one  often  would  not  cause 
Newton  to  converge.  Further,  when  convergence  of  Newton  did  ensue,  many  iterations 
were  usually  required  and,  of  course,  evaluating  B ^ Involves  a significant  amount 
of  computation.  Computational  results  (see  section  7)  show  a significant  Improvement 
with  the  new  logic.  Evidently,  once  the  radius  of  convergence  of  Newtons  method  (with 
B ^ evaluated  at  a = 0)  has  been  reached,  it  is  not  worthwhile  to  try  to  extend  it.  In 
fact,  it  had  proved  worthwhile  to  limit  motion  from  a * 0 to  those  o values  which 
require  only  a few  Newton  iterations  (fewer  than  6,  for  example)  for  convergence. 

This  is  accomplished  by  the  test  in  block  7,  page  2 of  the  DMINRG  flow  chart. 

Additional  advantages  of  this  new  logic  are  (1)  it  is  much  simpler  and  (2) 
only  one  evaluation  of  B ^ is  required  per  one  dimensional  search.  This  latter 
feature  is  especially  Important  for  large  problems,  where  evaluating  B ^ will  be  one 
of  the  major  computational  steps. 

Subroutine  REDOBJ  «ow  includes  logic  to  deal  with  basic  variables  which  violate 
their  bounds.  This  permitted  the  elimination  of  subroutine  BSCHNG  of  reference  [1], 
and  hence  shorted  the  code. 

The  logic  in  blocks  6 and  11  on  page  2 was  added  as  part  of  the  changes  required 
to  solve  problems  for  which  the  initial  solution  was  not  feasible  i.e.  a Phase  I to 
determine  an  initial  feasible  solution.  These  blocks  deal  with  constraints  which 
had  previously  violated  an  upper  (lower)  bound  and  now,  because  too  large  a step  has 
been  taken,  violate  their  lower  (upper)  bounds. 
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The  logic  on  page  A,  block  15,  which  stores  the  point  with  lowest  objective 
value,  is  also  new.  This  causes  each  one  dimensional  search  to  yield  a point  no 

worse  than  the  one  it  starts  with.  This  logic  was  added  because  the  last  point 
obtained  by  DMINRG  is  not  always  the  best,  especially  when  it  is  obtained  by 
interpolation.  This  may  be  due  to  the  Inacurracy  with  which  F(x)  is  obtained; 
interpolation  seems  to  produce  poorer  results  as  EPNEWT  is  increased.  Without  this 
logic,  function  values  obtained  at  the  end  of  successive  DMINRG  calls  sometimes 
increased,  especially  near  the  minimum  and  with  larger  values  of  EPNEWT.  One  consequen( 
of  this  was  that  many  iterations  could  take  place  without  the  objective  improving 
and  without  the  stop  criteria  in  subroutine  GRG  being  met.  The  new  logic  eliminates 
this  problem. 

Subroutine  CONSBS 

The  only  major  change  from  the  previous  version  is  in  the  way  degeneracy  is  handlei 
In  block  13,  p.  2 of  the  CONSBS  flow  chart,  the  N branch  is  taken  when  an  acceptable 
pivot  element  cannot  be  found  in  mode  2.  The  old  code  simply  printed  an  error 
message  and  stopped.  The  present  code  branches  to  page  A of  the  flow  chart,  which 
enters  a variable  at  its  bound  into  the  basis.  This  is  an  improvement,  but  the  logic 
is  still  not  s^'isfactory . Future  plans  for  dealing  with  degeneracy  are  discussed 


in  section  8. 
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7.  Computational  Results 

(a)  Comparison  with  Previous  GRG  Code 

The  changes  described  in  section  6 have  both  simplified  the  code  and  improved 
its  efficiency.  Table  1 below  gives  problem  characteristics  of  eight  test  problems, 
while  table  2 gives  comparative  results  using  the  new  GRG  logic  and  the  old  described 
in  reference  [1]  and  section  6.  While  the  number  of  one  dimensional  searches  has 
increased  slightly,  the  other  3 performance  measures  are  significantly  Improved. 


J 


NAME  OF 
PROBLEM 

NO.  OF 

VARIABLES 

NO.  OP 

EQUALITY 

CONSTRAINTS 

NO.  OF 

INEQUALITY 

CONSTRAINTS 

Kowalik- 

Osborne 

Quadratic 

4 

0 

3 

Colvil . 
Quadra L_c 

5 

0 

3 

Asaadi  Problem 
No.  4 

5 

3 

0 

Eight  Variable 
Spring 

8 

0 

23 

R.A.C.  Shell 
Dual 

15 

0 

5 

Seven  Variable 
Truss 

17 

10 

8 

RAC  Primal 

5 

0 

10 

GGP  Alklya- 
tion 

7 

0 

1 

14 

TABLE  1 - Characteristics  of  Test  Problems 
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PROBLEM 

1 

STATISTIC 

OLD 

NEW 

One  Dim- 
ensional 
Searches 

8 

8 

Newton 

Iterations 

135 

56 

Equivalent 
Function 
Calls  * 

283 

143 

BINV 

Computa- 

tions 

15 

8 

2 I 3 { 4 I 5 I t totals 

OLD  I NEW  OLD  I NEW  OLD  I NEW  OLD  INEW  |oLD|NW|  OLD  |NEW 


5 5 


8 9 519  no  87  48  543  412  87  78  1379  813 


11 

68 

75 

78 

1379 

813 

208 

2874 

1797 

11 

126 

85 

TABLE  2 - Results  of  New  Logic 

* Equivalent  Function  Calls  r.  Functions  Calls  + N*  Gr.'dient  Calls 

where  N = No.  of  variablis. 
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(b)  Comparison  with  Penalty  Methods. 


The  same  seven  test  problems  were  run  on  the  current  GRG  code  as  well  as  on 
the  interior  penalty  code  described  in  [11].  GRG  required  far  fewer  on  dimensional 
searches,  function  evaluations  and .gradient  evaluations.  While  some  of  this  reduction 
was  offset  by  the  requirement  of  matrix  inversion  and  solution  of  nonlinear  equations 
in  GRG,  it  was  noted  that  GRG  produced  more  accurate  solutions  for  most  of  the  problems. 
Differences  in  computation  times  (of  the  order  of  hundredths  of  a second)  could  not  be 
estimated  accurately  owing  to  the  masking  effects  of  mutlprogrammlng.  We  expect 
that  GRG  will  prove  superior  to  penalty  methods  for  large  problems  where  linear 
programming  technology  can  be  Implemented  very  effectively.  The  comparative 
performance  of  GRG  and  the  penalty  codes  is  shown  in  Table  3 


STATISTIC 


One  Dimensional 
Searches 


Function  Calls 


Gradient  Calls 


Equivalent 
Function  Calls 


Newton 

Iterations 


BINV 

Computations 


TOTALS  - PROBLEMS  1, 2,4-8  | 

PENALTY 

GRG 

495 

94 

3773 

1125 

539 

99 

REDUCTION  FACTOR 
PENALTY /GRG 


Table  3 - Comparison  of  GRG  and  Interior  Penalty  Codes, 
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8.  Future  Work 

a.  Degeneracy 

A feasible  point  is  degenerate  if  the  gradients  of  the  binding 
constraints  (including  bounds  on  the  variables)  are  dependent.  Hence,  a 
point  with  more  than  N constraints  binding  is  always  degenerate.  At  a 
degenerate  point,  any  basis  matrix  must  Include  at  least  one  column 
corresponding  to  a variable  at  one  of  its  bounds.  This  may  be  a slack 
variable  or  one  of  the  natural  variables.  As  one  attempts  to  move  away 
from  a degenerate  point,  it  may  be  that  one  or  more  of  the  basic  variables 
at  a bound  immediately  violates  that  bound.  Then  either  the  basis  must 
be  changed  or  a new  search  direction  selected.  We  have  encountered  3 
degenerate  points  in  the  dozen  or  so  problems  solved  thus  far.  The  current 
logic  for  dealing  with  degeneracy  is  unsatisfactory.  Work  in  the  coming 
year  will  aim  at  remedying  these  defects, 

b . Extrapolation  of  Basic  Variables 

Within  a given  one  dimensional  search,  each  basic  variable,  y^,  is 
a function  of  the  step  size  a,  y^  = y^(a) . At  a = 0 we  know  y^  (0)  and 
yj^  (0)  = V^,  where  is  the  i*"^  component  of  the  tangent  vector  V.  At 
a = B we  compute  y^^CB).  Given  this  information,  a quadratic  can  be  fit 
to  y^(0) , yJ^(O),  y^(B)  and  used  to  estimate  y^(C),  These  estimates  are 

used  as  initial  values  in  the  Newton  iteration  which  evaluates  y^(C). 

Then  a quadratic  can  be  fit  thru  y^^CO),  y^(B),  y^(C) . This  is  used  to 
estimate  y^  at  the  next  a value,  and  so  on.  At  each  stage,  we  have  a 
quadratic  approximation  to  the  curve  y(a)  on  the  constraint  surface.  These 
approximations  should  be  much  better  than  the  current  linear  approximation, 
which  always  uses  the  tangent  vector  V evaluated  at  a = 0.  Computational 
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experiments  to  measure  the  improvement  yielded  by  these  moving  quadratic 
extrapolations  will  be  carried  out  in  the  coming  year. 


c.  Conjugate  Gradient  Algorithm 

We  have  changed  from  the  PPP  to  the  BFS  variable  metric  algorithm, 
due  to  recent  work  by  Shanno  T9,  10,]  showing  it  to  be  superior.  However, 
all  variable  metric  algorithms  require  storage  and  updating  of  an  N by  N 
H matrix  which  is  not  sparse.  Hence  there  is  no  way  they  can  be  used  in 
a GRG  algorithm  that  is  to  solve  large  problems.  The  only  alternative 
algorithm  that  requires  no  H matrix  and  converges  finitely  on  quadratic 
functions  is  the  Conjugate  Gradient  (CG)  algorithm  [12].  We  plan  to 
Implement  a version  of  GRG  employing  CG  and  compare  its  performance  with 
that  of  the  BFS  procedure.  CG  should  permit  solution  of  problems  with  100 
constraints  and  two  to  three  hundred  variables  in  an  explicit  inverse  GRG 
code  which  does  not  exploit  sparsity. 

d.  Solution  of  Geometric  Programs 

Geometric  programs  (GP)  have  constraint  and  objective  functions 


of  the  form 


gk  (X)  = ^ H H (X) 

i="k-l 


where  the  terms  t^  are  given  by 

N 


t,(x)  = nxf  J 

j=i'’ 

The  exponents  a^j  are  arbitrary  real  numbers.  Such  problems  have  received  much 
attention  in  the  literature,  and  special  algorithms  based  on  the  dual  GP  have 
been  developed  to  solve  them.  These  algorithms  arc  not  completely  satjsfa-tory , 
however,  due  to  ill-conditioning 
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of  the  dual  objective  function  [14]  . We  Intend  to  solve  a number  of 
geometric  programs  using  GRG  and  compare  the  results  with  those  obtained 
by  other  researchers  with  idiom  we  are  In  correspondence.  Our  GRG  code 
will  be  specialized  to  the  GP  structure  by  adding  a "front  end"  module. 
This  will  permit  Input  of  the  GP  by  specifying  the  Indices  nj^,  constants 
Cj^,  and  exponent  matrix  a^^.  It  will  also  compute  derivatives  of  the 
functions  gj^  as  the  terms  are  evaluated.  This  will  speed  computation 
significantly,  and  will  yield  the  accuracy  of  analytic  derivatives  without 
any  user-supplied  PARSH  subroutine. 


e.  A GRG  Code  for  Large  Problems 

One  of  the  most  attractive  features  of  GRG  Is  that  It  Is  able  to 
accommodate  large  problems  (1000  or  more  constraints  and  many  more  variables) 
by  Incorporation  of  techniques  which  exploit  sparsity.  These  techniques 
are  extensively  used  In  modem  linear  programming  codes.  A long 
range  (2-3  years)  goal  of  this  research  Is  the  development  of  a large 
scale  GRG  code  and  Its  comparison  with  current  codes  for  solving  large 
NLP's,  e.g.  approximation  programming  codes  [13]  . 
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